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Abstract 

The cubic nonlinear Schrodinger equation with repulsive nonlinearity and elliptic 
function potential in two-dimensions models a repulsive dilute gas Bose-Einstein 
condensate in a lattice potential. A family of exact stationary solutions is presented 
and its stability is examined using analytical and numerical methods. All stable 
trivial-phase solutions are off-set from the zero level. Our results imply that a large 
number of condensed atoms is sufficient to form a stable, periodic condensate. 
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1 Introduction 



The cubic nonlinear Schrodinger equation with attractive or repulsive nonlin- 
earity and a potential is used as a mean-field model for the dynamics of a 
dilute-gas Bose Einstein condensate (BEC) [1,2]. In this case, the equation is 
often referred to as the Gross-Pitaevskii equation. These BECs are of interest 
in both the theoretical and experimental physics community: they are exam- 
ples of macroscopic quantum phenomena which display phase coherence [3-6] 
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and which have potential applications in such diverse areas as quantum logic 
[7], matter- wave diffraction [8], and matter- wave transport [9]. 



After rescaling of the physical variables, the governing equation is 




where A denotes the Laplacian operator, ip(x, t) is the macroscopic wave func- 
tion of the condensate in one, two or three dimensions, with x defined as x, 
(x,y) or (x,y,z) respectively. The real-valued function V(x) is an experimen- 
tally generated macroscopic potential. The parameter a determines whether 
(1) is repulsive (a = 1, defocusing nonlinearity) , or attractive, (a = — 1, 
focusing nonlinearity). Note that for BEC applications, both signs of a are 
relevant. 

Many BEC experiments use harmonic confinement, but recently there has been 
interest in confinement of repulsive BECs using standing light waves, resulting 
in a sinusoidal potential in the Nonlinear Schrodinger equation. So far, most 
interest has focussed on the quasi-one-dimensional regime, in which the longi- 
tudinal dimension of the condensate is far greater than the transverse dimen- 
sions, which are themselves of the order of the healing length of the condensate. 
In this case, the governing mean-field equation is (1), with a — 1, in one di- 
mension, and the experimentally generated potential V(x) = —Vq sin 2 (x). In 
[10,11], a number of exact solutions of this equation was constructed and their 
stability was investigated. In fact, a potential more general than a sinusoidal 
potential was considered: 



where sn(x, k) denotes the Jacobian elliptic sine function, with elliptic modulus 
k. In the limit k — 0, this potential reduces to a sinusoidal potential. It is 
important to realize that for values of k up to k = 0.9, the shape of the 
potential is virtually indistinguishable from a sinusoidal one. On the other 
hand, considering values of k close to 1, i.e., k > 0.999 gives periodic potentials 
with well-separated peaks, allowing the study of BEC dynamics in an entirely 
new regime. 

Currently, no experiments are being performed where a BEC is trapped in a 
higher-than-one-dimensional periodic potential. However, the interest in the 
applications mentioned above strongly suggests that these experiments may 
eventually take place. Already, there are theoretical investigations of BECs 
in multidimensional lattice potentials [12], suggesting the realization of such 
experiments. Although motivated by the developments in BECs, in this paper 
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Fig. 1. Various lattice potentials. For all figures, A\ = A 2 = B\ = B 2 = 1. For (a), 
ki = k 2 = 0, mi=m 2 = 1; for (b), h = 0.999, fc 2 = 0, mi = 2K(0.999)/vr, m 2 = 1. 
Finally, for (c) ki = k 2 = 0.999 and mi=m 2 = 2K(0.999)/vr. 

we consider (1) with repulsive nonlinearity in two dimensions in its own right. 
Thus we consider 

with if) = ip(x, y, t). As in one dimension [10,11], making the proper choices for 
the potential allows the construction of a large class of exact solutions. The 
family of potentials considered is 

V(x, y) = - (Aisn 2 (mi:r, h) + B^j (A 2 sn 2 (m 2 y , k 2 ) + B 2 ) + 

+mlklsn 2 (mix, k\) + m^k^sv? {m 2 y , k 2 ). (3) 

Here Ai, A 2 , B 1 , B 2 , mi and m 2 are real constants. The two elliptic moduli ki 
and k 2 are in the interval [0,1]. The first term is a straightforward general- 
ization of the one-dimensional potential (2), with an additive constant. The 
other terms are incorporated to facilitate the construction of exact solutions. 
Although this exact expression for the potential is necessary to allow the con- 
struction of exact solutions, it is the qualitative features, i.e., its periodicity 
and amplitude, that are most important. Numerical and analytical results 
throughout this paper demonstrate that the behavior of a solution in a lattice 
potential is largely independent of the quantitative features of the potential. 
Figure 1 displays the potential (3) for various values of its parameters. In the 
next section, a family of exact solutions to (1') is given and their linear sta- 
bility is discussed. Numerical results for various representative solutions are 
given in Section 3. The last section concludes the paper with a brief summary 
of the relevant mathematical and physical results. 
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2 Analytical Results 



A judicious choice of the potential allows for the cancellation of the nonlinear 
term in (1') so that exact solutions can be constructed. Of course, one can 
always find a suitable potential by solving (1') for V(x, y), given a certain 
ip(x,y,t). However, this generically results in time-dependent potentials and 
is hence not of interest. The kind of potentials we seek are two-dimensional 
generalizations of the one- dimensional potential (2). This one-dimensional po- 
tential is periodic. Its two-dimensional generalizations should be periodic in 
both its spatial variables. Furthermore, it is desirable that for a given potential 
more than one exact solution exists. 

A derivation of the exact solutions is not presented in this paper. Rather we 
state the solution formulas and discuss them. For the potential (3), a class of 
exact solutions is given by 

1>(x, y, t) = r 1 (x)r 2 (y)e i01 ^ +i0a ^- iut , (4) 

with 



/dz 
~A 2? Z'liP ' ( 5a ) 
| A\ smfmiz, k\) + B\ 

y 

/dz 
Amw-faua, ' (5b) 

and 



o 



\m\ (l + k\ + fc?^) + \m\ (l + k\ + kf^ , (6a) 
+ Bx) (A 1 + kfBA , (6b) 







UJ = 


\ ml ( 
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d = 


A 2 



+ B 2 ) [A 2 + k\B 2 ) . (6c) 

Choosing the parameters A±A 2 , Bi/A\, B 2 /A 2 , mi, m 2 , ki and k 2 determines 
the potential. Thus the solution class given by (5a-b) is a one-parameter family 
of solutions, with free parameter Ai/A 2 . Existence of these solutions requires 
rf > 0,cl> and r\ > 0, c 2 > 0. This imposes conditions on the parameters: 
Ai > 0, Bi > or Bi > 0, -A { < B { < -Ai/kf, where % = 1, 2. 

From (5a-b), it follows that ri(x) is periodic in x with period 2K(ki)/mi, 
where K(k) = f^ 2 dz/J 1 — k\ sin 2 (2), and r 2 (y) is periodic in y with period 
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2i^(/c 2 )/m 2 . In general, 9\{x + 2K{k\)/m 1 ) ^ 9i(x) + 2mr, for some integer 
n. Thus, the solution (4) is usually not periodic in x or in y. Imposing this 
periodicity requires a quantization of the phase 9i(x) and 9 2 (y) [11]. It is 
unclear if such a quantization is possible, since only one free parameter is 
available to satisfy two conditions. There are two special cases in which phase 
quantization is not a concern. The first case results in trivial-phase solutions, 
for which c\ = 0, c 2 = 0. The second case is the trigonometric limit, in which 
h x = 0, k 2 = 0. 

The solution has trivial phase in the i-direction if q is zero, where % — 1 (2) 
corresponds to the x (y)-direction. There are three possibilities: 

Bi = 0: n 
Bi = Ai : T'i 
Bi = Ai/ '. Ti 

where cia(miXi, fcj) is the Jacobian elliptic cosine function, and dn(mjXj, fcj) 
denotes the third Jacobian elliptic function. 

The solution is trigonometric in the i-direction if hi is zero. Then 



rf (x) = Ai sin 2 (xi) + Bi, tan(^) = y 1 + tan(m i a; i ), (8) 



'Ai sn(miX, hi), 



-Ai cn(rriiX, ki), 



-Ai/kf dnirriiX, ki), 



(7a) 
(7b) 
(7c) 



and phase quantization is satisfied. Notice that it is possible for the solution 
to have trivial phase in one direction and be trigonometric in the other. 

The main stability results from [10,11] can be generalized to two dimensions 
in a straightforward manner. The details of this will be presented elsewhere. 
Here a brief summary is given. Examining the linear stability of exact trivial- 
phase stationary solutions ip(x,y,t) = r(x,y) exp(-iut) of equation (1') gives 
rise to the eigenvalue problem 



( 












-:-) 










yu 2 ,) 



(9) 



where the linear Schrodinger operators L + and L_ are 



L + = - l - (d 2 x + dl) + 3r\x, y) + V(x, y) - 
L- = -\ {dl + dl) + r 2 (x, y) + V(x, y) - 



LO. 



(10a) 
(10b) 
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Reasoning similar to the arguments in [10,11] leads to the following conclu- 
sions: 

• If r(x, y) > 0, then r(x, y) is the ground state of L_ and ip(x, y, t) is linearly- 
stable. 

• If r(x,y) has any zeros, it is not the ground state of L_. If in addition L + 
is a positive operator, then ip(x,y,t) is unstable. 

• In all other cases, no conclusion is obtained from this stability analysis. 

Primarily, this analysis immediately establishes the linear stability of a solu- 
tion ijj(x,y,t) = r\(x)r 2 (y) exp(-iut), where both r\(x) and r 2 (y) are strictly 
positive. Thus, the solution with both the x- and y-dependence specified by 
(7c) is linearly stable. 



3 Numerical Results 

In this section, the result of numerically solving (1') with initial conditions 
chosen from the exact solutions given in the previous section are discussed. 
The numerical procedure uses a fourth-order Runge-Kutta method in time 
and a filtered pseudospectral method in space. For each experiment, a small 
amount of white noise was added to the initial conditions as a perturbation. 

Consider the solution 



ip(x,y,t) = -dn(m 1 x,k 1 )dn(m 2 y,k 2 )e lu)t . (11) 

where A\ < and A 2 < 0. As stated in the previous section, this solution 
is linearly stable. This is confirmed by the numerics presented in Fig. 2. The 
three different columns give from left to right the dynamics of the solution, 
the same using contour plots, and the arctan of the Fourier power spectrum 
of the solution. In the first two columns, the bottom figure shows the poten- 
tial V(x, y). The stability of this solution is reminiscent of the stability of a 
plane wave solution of the two-dimensional defocusing nonlinear Schrodinger 
equation [13]. 

Next, consider the trivial-phase solution 

i/j(x, y, t) = yf-A iy J-A 2 cn(m lX , fci)cn(m 2 y, k 2 )e~ iwt . (12) 

The numerics of Fig. 3 shows that this solution is unstable. For the equa- 
tion (1'), the L 2 -norm of the solution is conserved, hence so is the L 2 -norm 
of its Fourier spectrum. This is not reflected in Fig. 3, due to the arctan 
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Fig. 2. The stable solution from Eqn. (11). Here k\ = k 2 = 1/2, m\ = m 2 = 1 and 
Ai = A 2 = -1. 

transformation, which is used to diminish the range of the power spectrum. 
The onset of instability occurs between t — 15 and t — 30. A detailed look 
at this onset is given in Fig. 4. The middle column of this figure shows 
\\ip(x, y, t)\ 2 — \ip(x, y, t = 0)| 2 |. The instability begins to develop around t = 18. 
Before the instability occurs, \i/)(x, y, t)\ 2 is equal to its initial condition, up to 
effects due to the added noise and numerical round-off, as is expected, since 
ip(x,y,t) is an exact stationary solution of (1'). The solution 

ip(x, y, t) = yj ' A l A 2 sa(miX, k l )sn(m 2 y , k 2 )e~ iU}t . (13) 

is observed numerically to be unstable as well. Its dynamics is very similar to 
that of the solution (12). These results are consistent with the one-dimensional 
results obtained previously [10,11]. 
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Density: I ^(x^t) 1 2 Contour: l^(x,y,t)l 2 Fourier modes 




Fig. 3. The unstable solution from Eqn. (12). Here k\ = k 2 = 1/2, mi = m 2 = 1 
and Ai = A 2 = 1. 

Analytical results for the stability of the nontrivial phase solutions are difficult 
to obtain, and so far numerical investigations are the only means of examining 
these solutions. Since phase quantization is a serious complication for the 
numerics, it is convenient to work with solutions (8) that are trigonometric in 
both directions since phase quantization is automatically satisfied. Thus, we 
consider solutions of the form 



^(x,y,t) = ^(Aisin 2 (mia;) + (A 2 sm 2 (m 2 y) + B 2 )e ie ^ +ie ^- lw \ (14) 
with 



tan(0i(z)) = y 1 + -^tan(mix), tan(0 2 (y)) = \j 1 + -^tan(m 2 y). (15) 
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Fig. 4. A blow-up of the onset of instability occurring in Fig. 3. The middle figure 
shows \\ip(x,y,t)\ 2 - \ip(x,y,t = 0)| 2 | 

These solutions are stable or unstable depending on the offset parameters B\ 
and B 2 . As shown in Fig. 2 the off-set dn(mjXj, fcj)-type solution is stable 
whereas in Fig. 3 the unstable sn(rriiXi, hi) and cn(mjXj, fcj) solutions without 
offset are unstable. These trivial phase solutions suggest that offset is essential 
for stability. In Fig. 5, A 1 — A 2 — 1 and B 1 = B 2 = 0.5 and the onset of 
instability of the nontrivial phase solution is observed. Here the chosen offset 
(Bi — B 2 — 0.5) is insufficient to stabilize the dynamics. In contrast, Fig. 6 
has parameter values A 1 — A 2 — 1 and B 1 — B 2 — 1 and is stable. Thus 
with a higher offset, the nontrivial phase solution is stabilized. The boundary 
between the stable and unstable regions is difficult to calculate analytically and 
we rely on numerical simulations to determine the amount of offset required for 
stability. This behavior is again consistent with the one-dimensional case [11]. 
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Density: I ^(x^t) 1 2 Contour: l^x^t)! 2 Fourier modes 




Fig. 5. An unstable solution from Eqns. (14) and (15). Here k\ = k% = 1/2, 
mi = m2 = 1, Ai = A2 = 1, and B\ = B2 = 0.5. 

4 Summary and Conclusions 

We considered the repulsive nonlinear Schrddinger equation with an elliptic 
function potential in two dimensions. Periodic solutions of this equation were 
found and their stability was investigated both analytically and numerically. 
Using analytical results for trivial-phase solutions, we showed that solutions 
with sufficient offset are linearly stable. This is confirmed with extensive nu- 
merical simulations on the governing nonlinear equation. Likewise, nontrivial- 
phase solutions are stable if they are sufficiently off-set. 

Since our equations is a model for a Bose-Einstein condensate trapped in a 
lattice potential, our results imply that a large number of condensed atoms 
is sufficient to form a stable, periodic condensate. Physically, this implies sta- 
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Density: I ^(x^t) 1 2 Contour: I W(x,y,t)l 2 Fourier modes 




Fig. 6. A stable solution from Eqns. (14) and (15). Here k± = k 2 = 1/2, 
mi = 7U2 = 1, Ai = Ai = 1, and B\ = B 2 = 1. 

bility of states near the Thomas-Fermi limit. The results are consistent with 
the one- dimensional trapping of a BEC condensate in a standing light wave. 
To quantify this phenomena, we consider the k = limit and note that in the 
trigonometric limit k — > the number of particles n per potential well is given 
by n = (m^y,t)\ 2 ) = (r?(*)) (r 2 2 (y)) = (A 1 /2 + 5x)(A 2 /2 + B 2 ), where (•) 
denotes the spatial average. In the context of the BEC, and for a fixed atomic 
coupling strength, this means a large number of condensed atoms n per poten- 
tial well is sufficient to provide an offset on the order of the potential strength. 
This ensures stabilization of the condensate. Alternatively, a condensate with 
a large enough number of atoms can be interpreted as a developed condensate 
for which the nonlinearity acts as a stabilizing mechanism. 

Acknowledgments: The authors acknowledge useful conversations with J. C. 
Bronski, L. D. Carr, R. Carretero-Gonzalez, K. Promislow, and W. Reinhardt. 
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